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Abstract 

Oh We derive explicit lower and upper bounds for the probability generating functional 

of a stationary locally stable Gibbs point process, which can be applied to summary 
statistics like the F function. For pairwise interaction processes we obtain further esti- 
mates for the G and K functions, the intensity and higher order correlation functions. 

The proof of the main result is based on Stein's method for Poisson point process 
approximation. 
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^ ; 1 Introduction 



Gibbs processes are very popular point process models which are extensively used both in 
spatial statistics and in statistical physics, see e.g. j9, 13]. Especially the pairwise interaction 



processes allow a simple yet flexible modelling of point interactions. However, a major 
drawback of Gibbs processes is that in general there are no analytic formulas available for 
their intensities or higher order correlation functions. In a recent couple of articles 0, 0] 
Baddeley and Nair proposed an approximation method that is fast to compute and accurate 
as verified by Monte Carlo methods. There are however no theoretical results in this respect 
and hence no guarantees for accuracy in most concrete models nor quantifications of the 
approximation error. 
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The aim of the present paper is to derive rigorous lower and upper bounds for correlation 
functions and related quantities. These allow us to narrow down the true values quite 
precisely if the Gibbs process is not too far away from a Poisson process. Figure [TJ shows our 
bounds on the intensity for a two-dimensional Strauss process in dependence of its interaction 
parameter 7. The pluses are estimates of the true intensity obtained as averages over the 
numbers of points in [0, l] 2 of 10,000 Strauss processes simulated by dominated coupling from 
the past. The point processes were simulated on a larger window in order avoid noticeable 
edge effects. All simulations and numerical computations in this paper were performed in 
the R language [12| using the contributed package spat st at |4|. 
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Figure 1: Bounds on the intensities of two-dimensional Strauss processes with = 50, 
r = 0.05 and values of 7 ranging from to 1. The pluses are estimates of the intensities 
based on 10,000 simulations each. 

Our main result, Theorem 12 .![ more generally gives bounds on the probability generating 
functional of a Gibbs process. Let H be an arbitrary point process on M. d . The probability 
generating functional (p.g.fl.) is defined as 

* a (p)=E(j]^(y)) (1.1) 

yes 

for any measurable function g: JBL d — > [0, 1] for which 1 — g has bounded support, see e.g. js), 
p. 59] for details. 

Many statistics of point processes, such as the empty space function (F function), contain 
expectations as in (11.11) . For pairwise interaction processes the situation is even better. 
By the Georgii-Nguyen-Zessin equation (12.21) the nearest neighbour function (G function), 
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Ripley's K function and the correlation functions of all orders can be rewritten using the 
p.g.fl. 

The idea for proving Theorem 12.11 is to replace the Gibbs process S in (11. ip by a suitable 
Poisson process and bound the error using Stein's method. 

The rest of the paper is organised as follows. In Section [2] we introduce some notation 
and state the main result. In Section [3] we provide bounds on the intensity, and in Section H] 
bounds on other summary statistics are derived. Section contains the proof of the main 
result. 

2 Preliminaries and main result 

Let (9t,jV) denote the space of locally finite point measures on W 1 equipped with the a- 
algebra generated by the evaluation maps [0T 3 £ i-> G Z + ] for bounded Borel sets 

A G R d . A point process is just a Ot-valued random element. We assume the point processes 
to be simple, i.e. do not allow multi-points. Thus we can use set notation, i.e. x G £ means 
that the point x lies in the support of the measure £. 

In spatial statistics point processes are usually defined on a bounded window W C 
R d . Let 911 w denote the restriction of the 91 to W. A point process H on W is called a 
Gibbs process if it has a hereditary density u with respect to the distribution of the Poisson 
process with unit intensity. Hereditarity means that «(£) > implies u(rj) > for all 
subconfigurations rj C £. By hereditarity we can define the conditional intensity as 

A(x | = m(0 , (2.1) 

where 0/0 = 0. Roughly speaking, the conditional intensity is the infinitesimal probability 
that H has a point at x, given that S coincides with the configuration £ everywhere else. 
Furthermore A(- | •) uniquely characterises the distribution of H, since by ( 12. ip one can re- 
cursively recover an unnormalised density. It is well-known that the conditional intensity 
is the dx ® jSf (H)-almost everywhere unique product measurable function that satisfies the 
Georgii-Nguyen-Zessin equation 

El f h{x,E\{x}) S(dar) ) = / E(h(x, S)A(x | 5)) (2.2) 

for every measurable /i: W x — >• M + . 

So far A(- | •) is only a function on W x 9X|>v, but in many cases there exists a natural 
extension to the whole space, which we shall also denote by A(- | •). One way to generalise 
Gibbs processes to the whole space M d is then by the so-called integral characterisation. A 
point process H on IR d is a Gibbs process corresponding to the conditional intensity A(- | •) if 
it satisfies (2.2) with W replaced by R d for all measurable h: R d x 01 -> R + : sec p. 95]. 
or [ll| for a more rigorous presentation. Unlike in the case of a bounded domain, S may 
not be uniquely determined by ( 12. 2p . For the rest of this paper we will only deal with the 
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conditional intensity, i.e. if we say that a result holds for a Gibbs process with conditional 
intensity A(- | •), we mean that it holds for all processes corresponding to this conditional 
intensity. 

A Gibbs process E is said to be a pairwise interaction process if its conditional intensity 
is of the form 

X(x\0 = PU^(x,y) (2.3) 

for a constant (3 > and a symmetric interaction function (p. We denote the distribution of 
5 by PIP (/?,<£>). The process E is called inhibitory if (p < 1 and it is said to have a finite 
interaction range if 1 — 99 is compactly supported. E is stationary if <p(x,y) depends only 
on the difference x — y; we then write ip(x,y) = tp(x — y). If in addition the interaction 
function is rotation invariant, i.e. <p(x) = ip(\\x\\), then E is called isotropic. For conditions 
on </? ensuring the existence of E we refer the reader to [l3[ . 

If H is a general point process, its expectation measure or first order moment measure 
EE on R d is simpy given by (EE) (A) = E(E(A)) for every Borel set A C R d . For k e N the 
fc-t/i order factorial moment measure of H is the expectation measure of the factorial product 
measure 

SW = 6 (Xu-,x k ) 

pairwise different 

on (M d ) fc . Any moment measure is said to exist if it is locally finite. 

The intensity (function) X(x) of a Gibbs process 5 is the density of the first moment 
measure of E with respect to Lebesgue measure, provided the first moment measure exists. 
For a bounded A C M d , Equation ( |2T2|) yields 



ES(A) = / E(A(x|S)) tfa, 



hence the existence of the intensity and the form X(x) = E(X(x | 5)). For stationary processes 
the intensity is obviously constant and we just write A. For a stationary pairwise interaction 
process we get 

A = E(A(0|H))=/3E(jJv(y)) =/3*a(v). (2.4) 

j/gE 

In a similar manner it is possible to obtain the densities of the higher order factorial 



moment measures, the so-called correlation functions; see [10|, |8j. For a stationary process 
E ~ PIP (/J, y>) the /c-th correlation function is given by 



x k (x u . . . ,x k ) = /3 k (^ Yl vfa - a; i)) E (n v9 ^ ~ Xi ^> ' ' m( p(y ~ Xk ^) 

l<i<j<k y£5- 

= /3 fc ( J] ^(x l -x j )^^(--x l )---^(--x k )). (2.5) 



l<i<j<k 
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A frequently used function in spatial statistics is the pair correlation function which is defined 

as 

P(x,y) = w rw v (2-6) 
\{x)\(y) 

In the stationary isotropic case this simplifies to p(s) = \2(x,y)/\ 2 , where s = \\x — y\\. 

For our results we need a stability condition for the Gibbs processes. A Gibbs process S 
is called locally stable if there exists a non-negative function c* such that f w c*(x) dx < oo 
for all bounded domains W C W 1 and the conditional intensity satisfies 

\(x\Z)<c*(x), (2.7) 

for all fetft. 

For the rest of the paper we restrict ourselves to stationary point processes on JBL d and 
require c* to be a constant. The following is the key theorem for obtaining the results in 
Sections [3] and HI Its proof is the subject of Section [5j 

Theorem 2.1. Let H be a stationary locally stable Gibbs process with intensity A and local 
stability constant c* , and let g : M d — > [0, 1] be a function for which l—g has bounded support. 
Then 

1 - AG < E(n*(l/)) < 1 " - e " C * G )> (2-8) 
yes 

where G = L d 1 — g(x) dx. 



3 Bounds on the intensity 

For the intensity of a inhibitory pairwise interaction process we immediately obtain from 
Theorem 12.11 the following result. 

Theorem 3.1. Let 5 ~ PIP(/3,<£>) be inhibitory with finite interaction range. Then 

! <A<^-^, (3.1) 



1 + 0G ~ ~ 2 - e-^ 
where G = f Rd 1 — <p(x) dx. 

Proof. Recall from (12. 4p that A = 0^Y[yesf(y) anc ^ use c * = P- Theorem 12.11 then yields 
which can be rearranged as (13. ip . □ 



Remark 3.2. The lower bound of (13. ip can also be found in 13j, p. 96] with the restriction 
f3G < e _1 , whereas our inequality holds for all values of f3 and G. 
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Example 3.3. Let S be a Strauss process, i.e. 



if < r 

if ||a;|| > r 



for some parameters r > and < 7 < 1. Then G = (1 — r y)a, l ir d , where denotes 
the volume of the unit ball. Figure [T] shows that for a reasonable choice of the parameters 
(/3, r, 7) the bounds on A are quite good. The maximal relative error between the bounds 
and the simulated values is about 3.5%. 




Figure 2: Intensities of two-dimensional Strauss processes with (3 = 100, r = 0.05 and values 
of 7 ranging from to 1. The solid line is Xps, the dashed line is Xmf, and the grey area 
corresponds to the bounds in (13.11) . The pluses are estimates of the intensities based on 
10,000 simulations each. 



For a comparison of our bounds on the intensity to known approximations from the 
literature we concentrate on two methods. 

The first one is the Poisson-saddlepoint approximation proposed in 0]. The authors 
replaced in H2.4[) the Gibbs process H ~ PIP(/3, <p) by a Poisson process H\ PS with intensity 
Xps such that the following equality holds 



A P5 = EA(0|H Aps ) = /3e( H p(y)). (3.2) 
Solving this equation yields 



*™ - ™, (3.3) 
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where W is Lambert's W function, the inverse of x i— > xe x , and G = J Rd 1 — y?(x) drr as above. 

The second method is the mean-field approximation that was also described in 0] and is 
given by 

Am, = ™, (3.4) 

where T = — j Rd \og((p(x)) dx. Figure [2] shows the two approximations and our bounds from 
Inequality ( 13. ip for two-dimensional Strauss processes. 

In ji[ it is shown that under the conditions of Theorem 13.11 we have A > Xmf- The 
authors also conjectured, based on simulations for Strauss processes, that X PS is an upper 
bound for A. However, the next example indicates that this is not generally true. 
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Figure 3: The two processes in Example I3.4I restricted to the unit square. In the left panel 
is a realisation of the hard annulus process having 489 points, whereas Xps = 295.2. In the 
right panel is a hard core process realisation having 188 points. 



Example 3.4. Consider the process H ~ PIP(/3,<£>) with the interaction function 

if \\x\\ < r 
(f(x) = ■{ if r < ||x|| < R 
if \\x\\ > R 

for constants < r < R. We refer to this as a hard annulus process. It is a special case of 
a so-called multiscale Strauss process, see @, Ex. 6.2]. Let d = 2, (3 = 3000, r = 0.05 and 
R = V2r. Then 

Amf = 0, -A- = 122.1, Xp S = 295.2 and - & - m = 1500. 

1 + pG 2 — e~P G 
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An estimate of the intensity based on 300 simulations gave A = 493.8 > Xps- For comparison 
we also estimated the intensity of a Strauss hard core process (7 = 0) with the same (3 and 
G and obtained A = 193.3. Figure [3] shows that although the two processes have the same 
/3 and G, their realisations look quite different. All simulations were performed by long runs 
(10 7 steps) of Markov Chain Monte Carlo. 

We were not able to prove that A > Xps in this case, but bring forward the following 
heuristic argument for the observed phenomenon. The simulations shows that for large /3 the 
points tend to cluster on "island" of radius < r/2 which are separated by a distance > R. 
Since the points within each island do not interact, we expect the intensity to grow linearly 
in (3 for large 0. However Xps only grows logarithmically for large /3, so that at some point 
the intensity will overtake. 

Even if \p$ may not serve as a bound on A it remains useful as an approximation. 
Empirically its values stay relatively close to the simulated values, whereas the difference of 
our upper and lower bounds in (13. ip increases for large (3G. 

The following result in connection with Theorem 13.11 gives an upper bound on the error 
in Poisson saddlepoint approximation. 

Proposition 3.5. Under the conditions of Theorem \3.1\ we have 

< A P5 < -Z ZTjfi- ( 3 - 5 ) 



1 + /3G ~ " ~ 2-e-? G ' 
Proof. Since \p$ = W((3G)/G, it suffices to show the following two inequalities: 

< W(x) and W(x) < - X _ (3.6) 



1 + x ~ y ' v ' ~ 2 - e 

for all x > 0. The first one follows from x/(l + x) < log(l + x), see [l, Eq. 4.1.33], by 
transforming it to 



exp ( ) < x 

1 + x v l + x' ~ 

and applying the increasing function W on both sides. For the second inequality note that 



log(2 - e"") < — (3.7) 
2 — e x 

This holds because we have equality for x = and it is straightforward to see that the 
derivative of the left hand side is less than or equal to the derivative of the right hand side 
for all x > 0. A similar transformation as above and applying W on both sides again gives 
the second inequality in (13. 6p . □ 
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4 Summary statistics 



For a point process S the empty space function or F function is defined as the cumulative 
distribution function of the distance from the origin to the nearest point in H, i.e. 

F(t) = F(3y EE: \\y\\ < t) = 1 - P(S(B(0, t)) = 0) 

= 1 - E(J[ l{y i 1(0, t)}) = 1 - * s (l{- i B(0,t)}), 
yes 

where B(x, t) denotes the closed ball centred at x G JR d with radius t > 0. Thus for a locally 
stable process E with constant c* we obtain from Theorem 12.11 

— (1 - exp(-c*a d t d )) < Fit) < \a d t d . (4.1) 

c* 

Note that for a Poisson process with intensity A we may choose c* = A, in which case the 
lower bound in (14. ip is exact. A minor drawback of the bounds in (14. ip is that the intensity 
is in general not known and has to be estimated as well, e.g. by the methods of Section |3j 

The nearest neighbour function or G function is defined as the cumulative distribution 
function of the distance from a typical point of H (in the sense of the Palm distribution) to 
its nearest neighbour. For pairwise interaction processes H ~ PIP(/3, (p) the G function is 
computed in |S|, Sec. 5] as 

G(t) = 1 - ^E(JI t{y i mt)Mv)) = 1 " jMH- i B(0,f)MO)- (4.2) 
yes 

Thus if S is inhibitory and has finite interaction range, setting c* = in Theorem 12.11 yields 

2 - f - exp(-/3G t ) < G{t) < 1 - ^ + f3G t , (4.3) 

where G t = L d 1 — V 9 ( x )l{|l x ll > The left panel of Figure S] shows these bounds for the 

hard annulus process from Example 13.41 with parameters (3 = 70, r = 0.025 and R = 0.035. 
Let us furthermore assume that H is isotropic. Then the K function is defined as 

K{t) = a d d [ s d_1 p(s) ds, 
Jo 

where p is the pair correlation function. By (12. 6p . (12. 5 p and Theorem 12.11 we obtain bounds 
on p as 

- ^) < p(NI) < (j2 - x ( : ~ ex P(-^))) . ( 4 - 4 ) 

where = L d 1 — <p(y)<p(y — x) dy, and (in most cases numeric) integration of (I4.4p yields 
bounds on the K function. 
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Example 4.1. Let H be a Strauss process in two dimensions. Then 

— 7) 2 ( arccosf ) — 

" \ \2r J 27 



G x = 2nr 2 (l -7) - 2r\l 



\\x\\ 
~2r 



(4.5) 



see also [2J. Since we do not know the true intensity of the Strauss process, we plug in the 
bounds of ( 13. ip into (14.41) to obtain bounds on the K function. This procedure causes twice 
an error and therefore the estimates on K are good only for smaller values of (3G. The right 
panel of Figure H] shows these estimates for i3 = 40, r = 0.05 and 7 = 0. 
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Figure 4: Left: G function of the hard annulus process from Example 13.41 with parameters 
= 70, r = 0.025, R = 0.035. Right: K function of a Strauss hard core process with 
parameters (3 = 40, r = 0.05. In both panels the solid line is an estimate of the true function 
based on 1,000 simulations, the dashed line is the the true function for a Poisson process (in 
the case of the G function with the same intensity as the one obtained by simulation for the 
hard annulus process), and the grey area corresponds to the bounds computed in the text. 



5 Proof of Theorem 12.11 

The main strategy of the proof is to replace in (II. ip the process H by a Poisson process H, 
and then use Stein's method to bound the error 

E(n^))- E (n^))- (5-i) 

yGE j/GH 

In the context of Poisson process approximation, Stein's method for bounding expressions 
of the form |E/(S) — E/(H)| uniformly in / from a class of functions T was first introduced 
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in [5(. It has undergone numerous developments in the subsequent years, including a recent 



generalisation to Gibbs process approximation [14j. In what follows we briefly touch on the 



key points of Stein's method needed for our proof, referring to the nice expository chapter 



of [15j for details. 

Assume that g: ~R d — > [0,1] is such that A = supp(l — g) is compact. Consider the 
function /: 01 ->■ [0, 1], 

/(o=/«w=n»(y)- ( 5 - 2 ) 

For any £ e *3X\a let Z^ = {Z^{t)} t >o be a spatial immigration-death process on A with 
immigration rate v > and unit per-capita death rate, started in the configuration £. This 
is a pure-jump Markov process on 0^ that holds any state r] € OTl^ for an exponentially 
distributed time with mean + rj(A)), where \A\ denotes the Lebesgue measure of A; 

then a uniformly distributed point in A is added with probability + rj(A)), or a 

uniformly distributed point in r\ is deleted with probability r](A) / (u\A\ + 77(A)). 

The process has the Poisson process distribution on A with intensity v as its stationary 
distribution. If the process is started at the empty configuration 0, then Z$(t) is a Poisson 
process with intensity v{\ — e _< ) for every t > 0. Let {E^}^ be i.i.d. standard exponentially 
distributed random variables and introduce the death process D^(t) = ^ xe g t{E x > t}S x . 
Constructing Z$ and independently on the same probability space, Z^ can be represented 



as Zt(t) = Z 9 (t) + Dz(t) for every t > 0; see [15|, Thm. 3.5] 



Let H„ be a Poisson process with intensity v and define hf. 91 — )■ R as 

/>oo 

MO = M£U) = - / [E(/(%i(*))) - E(/(H,))] (5.3) 
>/ o 

By 15|, Lem. 5.2] /i/ is well-defined and satisfies 

E/(S) - E/(H„) =E [ [h f (S + S x ) - h f {E)] vdx + E [ [h f (E- 5 X ) - h f (E)] E{dx) 

J A J A 

= E / [hf(E + 4) - ME)] (f " A(a: | S)) dx 

J A 

= E / [^(S + 4) - M s )] (" - I S)) dx, (5.4) 



where we applied the Georgii-Nguyen-Zessin equation on W 1 to the function [(#,£) h-> 
lyi(x) (/i/(£) — + <y)l for obtaining the second equality. 
Equation ( 15. 4p is our starting point for further considerations. 

Proposition 5.1. Let E be a stationary Gibbs process with intensity A and conditional 
intensity A(- | •). Let g: M d — > [0, 1] be a function such that 1 — g has bounded support. Then 
for all v > 

E( J] <?(</)) = 1 - \ (1 - e- G ) + 1^), (5.5) 
ye- 



ll 



where 

I v {g) = z~ vG ^U e vGs (l-\[ (l-s(l-g(y)))) ds ■ Jjl-g(x))(X(x | S)-i/) dx ) . (5.6) 



Proof. It is well known that for the Poisson process H u we have 



E(n»H)^(-"/i-^ 



e"^; (5.7) 



see e.g. p, Eq. 9.4.17]. We then follow the main proof strategy laid out above in order to 
re-express 

E(n^))-^ G 

Starting from Equation (I5.4p . we may use the decomposition Z^ + $ x — + D$ x with 
independent and .D^ to see that for any £ G 

/■oo 

M£ + MO = / E/(^(t))-E/(^+i B (t))dt 

/■oo 

= / Ef(Z^(t))-Ef(Z $ (t)+D Sx (t))dt 
Jo 

/■oo 

= (1 -g{x)) / E/(Z € (t))P(D, a (t) ^ 0) dt 
Jo 



poo 

[l-g(x)) / Efface-* dt 
Jo 



1 - <?(*)) / E(f(Z 9 (t)) + (f(Zz(t)) - f(Z (t))) e-t dt. 



By further decomposing Z^ = Z$ + with independent Z0 and D^, we obtain 
E(/(Z^))-/(Z B (t)))=E( J] s(y))E( J] <?(y)-l) 

= exp ( - Kl - e-«)G) ( 11(1 - - <?(</))) - 1 

yet 

where for the first expectation we used (15 .7p and that Z$(t) is a Poisson process with intensity 
u{l — e~ f ); for the second expectation note that each point of £ survives independently up 
to time t with probability e _< . Thus in total 

= (1 - ^(x))e^ G r + e**^ (jj(l - e ~*(l - <?(y))) - l)]e~* rft 

= (1 - S^)) 1 ^^ + (1 - <?(*)K* G f e vGs (Hi 1 - - 9{y))) ~ l) ds 
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by the substitution s = e 1 . 

Plugging this into Equation (15. 4p and using EA(x | S) = A finally yields 

E( II 9(V)) - e- G = ^ (1 - e- G ) + 7„(<?). 

□ 

Proposition 5.2. Let 3 be a stationary locally stable Gibbs process with constant c* . Then 
for all < v < c* 

h(g) < L(g) <h(g), 

where 

IM = -^;( c *(! " e ^ G ) " Kl " ^ C * G )) < 0, (5.8) 
4(<?) = i(c*(l - e^ G ) - 1/(1 - e^ G )) > 0. (5.9) 
Furthermore I v {g) < /or a// z/ > c*. 

Proof. Since — z/ < A(x | S) — z/ < c* — v a.s., we get for v < c* the upper bound 



Ug) < (c* - v)e~» G G I e» Gs (l - E J] (l - «(1 - g(y)))) ds 

Jo yes 



(5.10) 



and a similar lower bound, where c* — v is replaced by —v. Because A = supp(l — g) is 
bounded, H can be replaced by in (15.101) . It is a known fact that every locally stable 
Gibbs process on a bounded domain can be obtained as a dependent random thinning of a 
Poisson process; see [3, Remark 3.4]. In particular, there exists a Poisson process H c . such 
that 5|^4 C H c * a.s. Since (1 — s(l — g(y)) < 1 for all s G [0, 1] and for all y G ~R d , we obtain 

1 - E J](l - a(l - <Ktf))) < 1 - E II (1 - '(I - g{y))) = 1 - e~ sc * G , 
yes y eH c , 

where the last equality follows by (15. 7p . Integrating and rearranging the terms yields the 
formulas (15. 8p and (15. 9p . If v > c*, I u (g) is obviously non-positive. □ 



Remainder of the proof of Theorem \2.1[ Propositions 15.11 and 15.21 yield the upper bounds 

_g(y)) < i- A 

yel 



E( II 9(V)) < 1 - £ (1 - + £ (1 " e- G ) - (1 - e-^) 



1 - p-^ G 
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for < v < c* and 



^(J[g{y))<l-\G- 



vG 

for v > c*. Since the function [x i— >■ (1 — exp(— is monotonically decreasing for x > 0, 
we obtain the minimal upper bound for v = c* , as 



yes 

For the lower bound recall the Weierstrass product inequality, which states 



n 

a, 

i=l i=l 



n(i-^)>i-E' 



for < ai,...,a n < 1. Then, noting that the products below contain only finitely many 
factors 7^ 1 by the boundedness of supp(l — g), we have 



E (n^^)= E (n( i -( 1 -^^)))! 

ye=. yes 

y(S 

= 1 — E / 1 - g{x) E(dx) 

= 1 — A / 1 - ^(a;) da; = 1 - AG 

by Campell's formula; see @, Section 9.5]. □ 

Remark 5.3. An alternative proof for the lower bound can be obtained by using Proposi- 
tions 15.11 and 15.21 in the analogous way as for the upper bound. 



References 

[1] Milton Abramowitz and Irene A. Stegun. Handbook of mathematical functions with 
formulas, graphs, and mathematical tables, volume 55 of National Bureau of Standards 
Applied Mathematics Series. 1964. 

[2] Adrian Baddeley and Gopalan Nair. Approximating the moments of a spatial point 
process. Stat, 1:18-30, 2012. 

[3] Adrian Baddeley and Gopalan Nair. Fast approximation of the intensity of Gibbs point 
processes. Electron. J. Statist., 6:1155-1169, 2012. 



14 



Adrian Baddeley and Rolf Turner. Spatstat: an R package for analyzing spatial point 
patterns. Journal of Statistical Software, 12(6):l-42, 2005. 

Andrew D. Barbour and Timothy C. Brown. Stein's method and point process approx- 
imation. Stochastic Process. AppL, 43(1):9-31, 1992. 

Daryl J. Daley and David Vere- Jones. An introduction to the theory of point processes. 
Vol. II. Probability and its Applications (New York). Springer, New York, second 
edition, 2008. 

Wilfrid S. Kendall and Jesper M0ller. Perfect simulation using dominating processes on 
ordered spaces, with application to locally stable point processes. Adv. in Appl. Probab., 
32(3):844-865, 2000. 

Shigeru Mase. Mean characteristics of Gibbsian point processes. Ann. Inst. Statist. 
Math., 42(2):203-220, 1990. 

Jesper M0ller and Rasmus P. Waagepetersen. Statistical inference and simulation for 
spatial point processes, volume 100 of Monographs on Statistics and Applied Probability. 
Chapman & Hall/CRC, Boca Raton, FL, 2004. 

Jesper M0ller and Rasmus P. Waagepetersen. Modern statistics for spatial point pro- 
cesses. Scand. J. Statist., 34(4):643-684, 2007. 

Xuan-Xanh Nguyen and Hans Zessin. Integral and differential characterizations of the 
Gibbs process. Math. Nachr., 88:105-115, 1979. 

R Development Core Team. R: A Language and Environment for Statistical Computing. 
R Foundation for Statistical Computing, Vienna, Austria, 2012. 

David Ruelle. Statistical mechanics: Rigorous results. W. A. Benjamin, Inc., New 
York-Amsterdam, 1969. 

Dominic Schuhmacher and Kaspar Stucki. On bounds for Gibbs point process approx- 
imation. Preprint, 2012. Available at http://arxiv.org/abs/1207.3096, 

Aihua Xia. Stein's method and Poisson process approximation. In An introduction to 
Stein's method, volume 4 of Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap., pages 
115-181. Singapore Univ. Press, Singapore, 2005. 



15 



